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We study the liquid-vapor phase behaviour of a polydisperse fluid using grand canonical simula- 
tions and moment free energy calculations. The strongly nonlinear variation of the fractional volume 
of liquid across the coexistence region prevents naive extrapolation to detect the cloud point. We 
describe a finite-size scaling method which nevertheless permits accurate determination of cloud 
points and spinodals from simulations of a single system size. By varying a particle size cutoff we 
find that the cloud point density is highly sensitive to the presence of rare large particles; this could 
affect the reproducibility of experimentally measured phase behavior in colloids and polymers. 
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Many complex fluids such as colloidal dispersions, liq- 
uid crystals and polymer solutions are inherently polydis- 
perse in character: their constituent particles have an es- 
sentially continuous range of size, shape or charge. Poly- 
dispersity is of significant practical importance because 
it can affect material properties in applications ranging 
from coating technologies and foodstuffs to polymer pro- 
cessing [1]. However, our understanding of the fundamen- 
tal properties of polydisperse fluids remains very limited 
compared with what we know about their monodisperse 
counterparts. The challenge arises because a polydis- 
perse fluid is effectively a mixture of an infinite number 
of particle species. Labelling each by the value of its poly- 
disperse attribute cr, the state of the system (or any of 
its phases) has to be described by a density distribution 
p{cr), with p{a)da the number density of particles in the 
range a . . .a + da. The most common experimental situ- 
ation is that in which the form of the overall or "parent" 
distribution p'^{a) is fixed by the synthesis of the fluid, 
and only its scale can vary depending on the proportion 
of the sample volume occupied by solvent. One can then 
write ^"((7) = n°/°((T) where /°(cr) is the normalized 
parent shape function and = N/V the overall parti- 
cle number density. Varying at a given temperature 
corresponds to scanning a "dilution line" of the system. 

A central issue in the physics of polydisperse fluids is 
the nature of their phase behaviour: in order to process a 
polydisperse fluid one needs to know under which condi- 
tions it will demix and what phases will result. However, 
the phase behaviour of polydisperse systems can be con- 
siderably richer than that of monodisperse systems, due 
to the occurrence of fractionation [2-4]: at coexistence, 
particles of each a may partition themselves unevenly 
between two or more "daughter" phases as long as-due 
to particle conservation-the overall density distribution 
p^{(j) of the parent phase is maintained. As a conse- 
quence, the conventional fluid-fluid binodal of a monodis- 
perse system splits into a cloud curve marking the onset 



of coexistence, and a shadow curve giving the density of 
the incipient phase; the critical point appears at the in- 
tersection of these curves rather than at the maximum of 
either [5]. 

In this letter we describe a joint simulation and theo- 
retical study of a model polydisperse Lennard- Jones fluid 
in which the size of the particles influences not only the 
length-scale but also the strength of the interparticle 
potentials, as defined in (2) below. For the case of size- 
independent interaction strengths, the critical point lies 
very close to the maximum of the cloud curve [6] , whereas 
for the present model we find that it is substantially be- 
low (see Fig. 3), as is observed in many experiments on 
complex fluids (see e.g. [7]) and simphfied theoretical cal- 
culations [8]. At the critical temperature, Tc, there then 
exists a finite density range where phase separation oc- 
curs on the dilution line. Most results shown below will 
be at this temperature; note that we will be interested 
mainly in the low-density part of the coexistence region 
rather than the critical effects at the other end, using Tc 
merely as a convenient temperature scale. 

The simulations were performed within the grand 
canonical ensemble (GCE). This is particularly useful 
for polydisperse systems, where it permits sampling of 
many different realizations of the particle size distribu- 
tion while catering naturally for fractionation effects. 
Operationally, we ensure that the ensemble averaged den- 
sity distribution always equals the desired parent distri- 
bution p^{(t) by controlling an imposed chemical poten- 
tial distribution p{a). A combination of novel and exist- 
ing techniques [9] are required to tune p{a) such as to 
track the dilution line, i.e. to vary the parent density n'^ 
but not its shape /°(cr). 

In the GCE simulations, the number density n is a fluc- 
tuating variable with average equal to vP. Its distribution 
p{n), shown in Fig. 1 for a range of parent densities n'^ 
a.t T — Tc, is a key observable. In the coexistence region 
it has two distinct peaks, which we sample using multi- 
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FIG. 1: Number density distribution p(n) at T = Tc for 
parent densities n° as indicated and for particle size cutoff 
ffc = 1-4. The system size is L = 15ff. See text after eq. (2)) 
for definitions of ctc and a. (a) Linear and (b) log scale. Inset: 
Liquid fractional volume v\ versus n°, for = 1-4, 1.6, 1.8. 



canonical preweighting [10]. The weight under the low 
and high density peaks corresponds respectively to the 
fractional volumes Wg and v\ that would be occupied by 
gas and liquid in the corresponding canonical ensemble. 
As expected, the peaks separate and the valley between 
them deepens as wc move away from the critical point by 
decreasing n". Concomitant with this is a gradual trans- 
fer of weight from the liquid to the gas peak. Finally the 
liquid peak disappears, at exponentially small values of 
v\ visible only on a log scale (Fig. 1(b)). 

The observed variation of p{n) raises the question of 
how to detect the cloud point n[!i, defined as the lowest 
parent density iiP where stable phase coexistence occurs. 
In a monodispcrsc system this is straightforward because 
the cloud point also gives the density of the gas phase, 
which remains constant throughout the coexistence re- 
gion. One then simply detects the point where the gas 
and liquid peaks have the same weight, i.e. r = v\/vg = 1, 
and measures the gas density there. (The criterion r = 1 
has the added advantage of leading to only exponentially 
small finite-size corrections to the value of /x at coexis- 
tence [11]). However, this method fails in a polydisperse 
system because fractionation causes the densities and size 
distributions of the coexisting phases to vary with [5] . 
One could attempt to locate the cloud point instead by 
extrapolating in to the point where v\ ^ [6]. But 
in our system the dependence of v\ on vP is so strongly 
nonlinear — another effect of fractionation, see inset of 
Fig. 1 — that the resulting cloud point estimates would 
have very large error bars. Indeed, on a linear plot of v\ 



versus n° as shown in Fig. 1(a) the effects of the particle 
size cutoff (7c which our more careful analysis will reveal 
(see Fig. 3 below) are essentially invisible. 

To make progress, wc analyse the finite-size scaling of 
p{n). As the linear system size L grows at fixed and 
T, the peaks in p{n) will narrow around the densities 
of gas and liquid, respectively, and the size distributions 
averaged over configurations from each peak will tend 
to those in the coexisting phases. The ratio r = 
is determined by the difference in the grand potential; 
this is directly related to the pressure P so that r = 
exp(/3i'^AP) for large L where jS = l/ksT and AP = 
Pi — Pg. The criterion for stable coexistence at given 
fixed nP is that r must have a finite value as L ^ oo; 
the pressure difference then has to scale as AP ~ L"*^ 
except in the special case r = 1 (see above). 




0.085 



0.075 0.08 0.085 0.09 0.095 

n« 

FIG. 2: Ratio r of liquid to gas fractional volumes on approach 
to the cloud point at T = Tc for ac = 1.4. The inset shows 
the (negative, scaled) second derivative of Inr w.r.t. n°. The 
peak position gives an estimate of the cloud point density. 
Squares indicate the scaled master curve (1). 

For finite L, mctastable coexistence can still be ob- 
served in the density region n° < n°i where AP = 0(1), 
but here r will be exponentially small. Fig. 2 shows 
this effect clearly: r is independent of L for sufficiently 
large nP, but the curves depart rapidly from each other 
(note the log-scale) for smaller nP. The cloud point sep- 
arates the two regimes, permitting the estimate nj^j w 
0.0825±0.0005 for the parameters shown in the figure. To 
derive a method which can estimate n°i even from data 
for only a single system size L, we use the fact that AP 
is 0(1) and scales linearly with n° — n°[ to leading order 
near the cloud point, and hence Inr ~ L'^{'nP — n^i). This 
applies for rP < n[![, while above one has Inr = 0(1). 
Thus the derivative (9/9n°)lnr should drop from an 
0{L'^) plateau to 0(1) around n° = n^y In the second 
derivative —{d/dn^ Y I'l'' this drop will manifest itself as 
a peak. A smooth derivative can be extracted from sim- 
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ulation data using histogram reweighting, and the peak 
position then serves as an estimate for n^j. This is shown 
in the inset of Fig. 2, and gives n^j ~ 0.0823 from the 
largest L, consistent with our earlier estimate derived 
from comparing data for different L. 

The above arguments can be formalized using the re- 
sults of [11], which pertain to the monodispcrsc case but 
which we have generalized to polydisperse systems [12]. 
We find that for large L the second-derivative plot ap- 
proaches a universal master curve 
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parameterized by z. The scaled parent density is defined 
here as n° = aL''-{inP — J^ci) + ln(6L''n[![) with a and h 
system-specific dimensionless scale factors. This scaling 
implies that the cloud point estimate from the peak posi- 
tion has finite-size corrections of order L^'^ In L, while the 
peak width and height scale as L^'^ and L^'', respectively. 
Our data are consistent with the width and height scal- 
ing and with the dominant L^'^ dopendonco of the peak 
shifts [12]. The master curve (1), appropriately scaled, 
is overlaid onto the largest-L data in Fig. 2 (inset) and 
shows excellent agrcicment. 

Fig. 1 shows that the metastable liquid peak in p{n) 
persists until well below the cloud point n^y The point 
at which it disappears marks the so-called mean field 
spinodal [13] where the liquid is unstable to small den- 
sity fluctuations. The parent density n^p where this oc- 
curs should tend to an L-indcpendcnt value as L grows 
large [13], and our data (not shown) are consistent with 
this. Spinodals in monodisperse systems are convention- 
ally characterized by the density of the phase that be- 
comes unstable, which is located inside the region where 
phase separation occurs. Here we use instead the den- 
sity rigp of the coexisting stable phase, which is outside 
this region. This is a more meaningful representation in 
the polydisperse context since only the stable (majority) 
phase has the parental size distribution, while that of the 
metastable (minority) phase is determined indirectly via 
chemical potential equality. 

Equipped with a systematic method for determining 
cloud points, we now consider the overall phase diagram 
of our system, the interparticle potential of which was 
assigned the Lennard-Jones form: 



(2) 



with 



cTiO-j, (^ij = (o'i + crj)/2 and rj. 



> 2.5(Tij and no 



The potential was truncated for 
tail corrections were applied. The diameters a are 
drawn from a (parental) Schulz distribution f^{a) oc 
a" exp [(z + 1)(t/(t], with a mean diameter a which sets 
our unit length scale. We chose z = 50, corresponding 



to a moderate degree of polydispersity: the standard de- 
viation of particle sizes is (5 = l/y/z + 1 14% of the 
mean. The distribution /°(cr) was limited to within the 
range 0.5 < a < ac- The upper cutoff ctc serves to pre- 
vent the appearance of arbitrarily large particles in the 
simulation, but would also be expected in experiment be- 
cause in the chemical synthesis of colloid particles, time 
or solute limits restrict the maximum particle size [14]. 

We complement the simulations with theoretical phase 
behaviour calculations, following closely our study of the 
purely size-polydisperse case [6] . An accurate expression 
for the excess free energy of a polydisperse hard sphere 
fluid accounts for the repulsive interactions. To this is 
added a van der Waals term which represents the attrac- 
tive part of the Uij . It scales as 



/ 



da da' p{a)p{a') {aa'){a + a'f 



(3) 



where the factors aa' and {a + a')^ arise, respectively, 
from the size dependence of the interaction amplitude 
eij and the interaction range Cy . Multiplying out 
gives an expression involving only the moment densities 
/ da p{a)a^ with i = 1 ... 4. Since the repulsive part of 
the excess free energy has a similar moment structure, 
the moment free energy (MFE) method [15] can be used 
for accurate numerical prediction of phase behaviour [6] . 
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FIG. 3: Comparison of cloud curves for CTc = 1.4 and CTc = 1.6. 
The critical points for tJc = 1.6 (x) and Uc = 1.4 (+) are 
marked. Also shown is the spinodal (limit of met ast ability) 
for (7c = 1.6. The inset displays the variation of the gas cloud 
point density n°i at T = Tc as a function of cTc, as obtained 
from GCE simulations (open symbols) and MFE theory (filled 
symbols). 

Fig. 3 shows cloud curves for upper size cutoff (Tc = 1.4 
and 1.6 as obtained from the GCE simulations. A strong 
(Tc-dependence is seen even though both values of ac are 
far in the tail of the parent distribution. This is at- 
tributable to very strong fractionation effects (Fig. 4): 
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FIG. 4: Size distributions /(cr) in the liquid shadow phase 
distributions a,t T — Tc for Uc = 1.4, 1.6, 1.8, together with 
the parent distribution /"(ct). Inset: MFE theory prediction 
for larger cutoff Uc = 3; note the second peak in the shadow 
distribution. 



despite particle sizes around CTc being very rare in the 
parent, they occur in significant concentration in the 
shadow liquid. Physically, since large particles interact 
more strongly, their presence leads to a substantial free 
energy gain at the shorter interparticle separations of the 
liquid. 

One is led to enquire whether the gas phase cloud point 
density would eventually tend to a nonzero limit as CTc 
is increased. The inset of Fig. 3 shows the simulation 
results and theoretical predictions. The former exhibit a 
further strong decrease of n^j by ^ 25% between ac = 
1.6 and 1.8; the latter suggest that this trend continues 
and that the cloud point density tends to zero for large 
(Tc. Such an unusual effect has previously been seen in 
theoretical studies of polydisperse hard rods with wide 
length distributions [16] and is also predicted to occur in 
solid-solid phase separation of polydisperse hard spheres 
[17], though only for large ctc and distributions with fatter 
than exponential tails. Here the decrease of n|?[ is clear 
even for CTc of 0(1), i-e. of the same order as ct, and 
scaling estimates suggest that cutoff effects occur for any 
size distribution with tails heavier than a Gaussian [12]. 

The physical origin of the decrease of n°i to zero is 
the appearance (for large ctc) of a second peak in the 
shadow phase size distribution near (Tc (Fig. 4, inset). 
As with the hard rods, we expect this second peak to 
eventually dominate as ac increases so that the shadow 
phase consists of ever more strongly interacting particles 
whose sizes are concentrated near ac- We speculate that 
as a consequence, there exists some cutoff for which the 
shadow phase liquid freezes into a quasi-monodispersc 
crystal phase. Indeed our simulations provide evidence 
for this scenario: for the large cutoff CTc = 2.8 the liquid 
spontaneously freezes to an f.c.c. crystal structure [12]. 
Although we observe this only for small n° values close 



to the spinodal point it is conceivable that, for ac values 
larger than those presently accessible to simulations, the 
freezing might occur from the stable liquid phase. 

Finally, with regard to the cloud curves as a whole 
(Fig. 3), we note that significant cutoff-dependent shifts 
occur only for densities below the critical density. This is 
consistent with our interpretation above: for higher den- 
sities, the shadow phase is a gas of lower density than 
the parent. In this, the concentration of large particles 
is suppressed and that of small particles negligibly en- 
hanced because of their weak interactions. The shadow 
size distributions are therefore concentrated well within 
the range 0.5 ... (Tc (data not shown) so that no cutoff 
dependence arises. 

In summary, the task of accurately locating cloud 
points of polydisperse fluids via simulation is severely 
complicated by fractionation effects. We have presented 
a generally applicable finite-size scaling method which 
addresses this problem. Application to a model poly- 
disperse fluid reveals the cloud curve to be highly sen- 
sitive to the presence of very rare large particles. Such 
effects imply that in experiments on polydisperse fluids 
(see e.g. [4]) it may be important to monitor and con- 
trol carefully the tails of the size (or charge, etc) distri- 
bution. Otherwise undetected differences could lead to 
large sample-to-sample fluctuations in the observed phase 
behaviour. 
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